Resummation of projectile-target multiple scatterings and parton saturation 
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In the framework of a toy model which possesses the main features of QCD in the high energy 
limit, we conduct a numerical study of scattering amplitudes constructed from parton splittings 
and projectile-target multiple interactions, in a way that unitarizes the amplitudes without however 
explicit saturation in the wavefunction of the incoming states. This calculation is performed in 
two different ways. One of these formulations, the closest to field theory, involves the numerical 
resummation of a factorially divergent series, for which we develop appropriate numerical tools. 
We accurately compare the properties of the resulting amplitudes with what would be expected if 
saturation were explicitly included in the evolution of the states. We observe that the amplitudes 
have similar properties in a small but finite range of rapidity in the beginning of the evolution, as 
expected. Some of the features of reaction-diffusion processes are already present in that range, 
even when saturation is left out of the model. 



I. INTRODUCTION 

Saturation of parton densities is a phenomenon that is 
expected to occur in the scattering of very fast hadrons 
P, 0]. It is the statement that the number of partons 
per unit of transverse phase space in the Fock states of 
the incoming hadrons does not grow exponentially for- 
ever with energy (or rapidity), as would come out of a 
naive solution of linear evolution equations: The growth 
gets softer at very high energies, in such a way that the 
unitarity of the probabilities of scattering be preserved. 

However, it seems very difficult to get saturation from 
a field-theoretical calculation. So far, the problem has 
not even been formulated properly in QCD, for it is al- 
ready very challenging to identify the graphs that should 
be taken into account. Equations such as the Balitsky 
equation [3[ (which is in fact a hierarchy of equations) 
have been written down, but they do not exhibit satura- 
tion in an obvious way (if at all), and they are anyway 
extremely difficult to solve. 

Despite these fundamental difficulties, new quantita- 
tive results could be obtained for QCD amplitudes in 
the saturation regime over the last few years [j, Q . In 
short, they rely on an analogy with reaction-diffusion 
processes [5| but not on the computation of definite Feyn- 
man graphs. The obtained results do not depend on the 
way how partons saturate, but they seem to require that 
there be such saturation phenomena. So saturation was 
assumed rather than found in these approaches. The 
matching of this statistical picture with field-theory cal- 
culations has been attempted Q , and has led to the state- 
ment that the Balitsky- Jalilian Marian-Iancu-McLerran- 
Weigert-Leonidov-Kovner (B-JIMWLK) equations 0, 0, 
U were incomplete, and that they had to be supple- 
mented with new terms. The problem was identified as 
follows: The Balitsky equations only contain splittings 
of partons, while mergings are needed in order to achieve 
saturation. The hope that merging rates could be ob- 
tained by boosting splitting vertices was turned down 
by the finding that this procedure would lead to nega- 
tive probabilities [JJ, [T3| , which is at least inconvenient in 
practice fllj . if not completely meaningless. 



In this paper, we go back to the original formulation of 
saturation by Mueller in the context of the color dipole 
model 0, [111, [3] , which was in fact based on the assump- 
tion that saturation of the parton densities is equiva- 
lent to unitarization of the scattering amplitudes through 
multiple exchanges between linearly-evolved Fock states 
of the incoming particles (if the rapidity is not too high) . 
Analytical calculations have been achieved within simi- 
lar approximations (see Refs. [HI, [IB] for recent progress), 
but the result looks always complicated in QCD and thus 
difficult to play with and to interpret. Numerical studies 
were conducted by Salam [13, EM 13 1 but at that time, 
there was no good theoretical understanding of the prop- 
erties that scattering amplitudes should exhibit at high 
energy. In the light of our present knowledge of satura- 
tion, that enables us to characterize saturation by ana- 
lyzing the properties of some traveling waves, we evaluate 
numerically, on toy models, how good this procedure is 
when only splittings and multiple exchanges are allowed. 
Another important highlight of the present work is that 
we are able to resum numerically the asymptotic series 
that can be constructed out of an expansion in the num- 
ber of rescatterings, and which has the structure of the 
proper field-theoretical series of successive Pomeron ex- 
changes. 

Our study assumes the following standard picture of 
scattering in the QCD parton model (see Fig. []}: An 
asymptotic hadron, made of valence partons, evolves into 
a set of quarks and gluons spread in impact parameter 
space, when its rapidity is increased. The rules of split- 
ting (and recombination when nonlinear effects are taken 
into account in the evolution) of the partons are fixed by 
the QCD Lagrangian. Two such hadrons interact with 
each other by exchanging, with probability of the order 
of a 2 (a is the strong coupling constant), one or several 
gluons between the partons of similar sizes and impact 
parameters that are present in the wavefunctions of the 
two hadrons at the time of the interaction. 

The outline goes as follows. In the next section 
(Sec. |TT|, we introduce a model for partonic state evolu- 
tion under splittings only. Sec. Mil is devoted to the for- 
mulation of the unitary scattering amplitudes built from 
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FIG. 1: Sketch of the scattering of the two evolved hadrons 
through splittings only (left) and through splittings supple- 
mented by nonlinear effects (here in the form of internal 
rescatterings; rightmost figure.) The vertical thick lines would 
correspond to gluon exchanges in QCD. Each of them comes 
with an extra power of a 2 which is not enhanced by corre- 
sponding powers of the rapidity. 



these states, that we compare to cases in which satura- 
tion is included explicitly in the evolution. The numerical 
study of the many variants of the toy model, in different 
frames and within different unitarization schemes, is the 
object of Sec. |TV] In Sec. |V] we present an alternative 
calculation, based on the numerical resummation of the 
asymptotic series in the number of Pomcrons that are ex- 
changed. We state our conclusions and some prospects 
in the last section. 



II. TOY MODEL FOR LINEAR PARTON 
EVOLUTION 

A. Construction of the model 

Focussing on one particular impact parameter, one 
may reduce the QCD problem to one transverse dimen- 
sion only. It is very important to keep the dynamics in 
the transverse space if one wants the model to be repre- 
sentative of some of the physics of QCD. However, we do 
not aim at incorporating all aspects of QCD: We need a 
simple model, that is easy to formulate in a way which 
can be implemented in the form of a Monte-Carlo event 
generator. Let us assume that there be only one type 
of object in the theory, which could be gluons or, more 
accurately, color dipoles, and that the latter may be fully 
characterized by a single "space" variable, which repre- 
sents for example their size in the transverse plane. There 
is in addition an evolution variable: the rapidity. So far, 
we have in fact just idealized the color dipole model 

To further simplify the model, we discretize it both in 
space and in rapidity. The evolution rule is the following: 
When the rapidity y is increased by one unit, a particle at 
position i on the lattice may be replaced, with probability 



5, by two offspring at respective positions i+j and i — j. 
For the distribution of j, we choose: 



Proba(j) 



1 



1 



(1) 



This rule obviously leads to an exponential increase of 
the number of objects on each site, which can eventu- 
ally break the unitarity constraints. Nevertheless, we do 
not a priori specify a saturation mechanism at this level: 
Several options of implementing unitarity of the scatter- 
ing amplitudes will be examined in the next section, and 
saturation, in the sense that the number of particles is 
prevented from exhibiting an exponential growth forever, 
will only be one of them. 

Note that the choice of discrete transverse space and 
especially discrete y is not very natural for a model that 
is meant to mimic QCD, in particular since discretizing 
y obviously breaks boost invariance (which makes sure 
that the rapidity evolution can be shared arbitrarily be- 
tween the two incoming hadrons without affecting the 
observables) . However, this choice is dictated by techni- 
cal reasons: We will need to be able to generate millions 
of events within some reasonable computer time. 

Let us define the generating function for the probabil- 
ity of the different configurations as 



Z(y,{xk}) 




P(y,{nk}), (2) 



where P{y, {p-kY) is the probability of having 
, rife, • • • particles on sites , k, ■ ■ ■ 

at rapidity y. This function contains all information 
about the statistics of the particle numbers on all sites. 
For example, the correlator of the number of particles 
on sites 1 and 2 is obtained by taking two derivatives: 



dZ 



dx\dx2 



(711712) 



(3) 



The generating function Z obeys the evolution equation 
Z(y + l,{x k }) = ±Z(y,{x k }) 

+ \( l - ~) E e ~' z (y> fa-AWn { W)- (4) 



' 3=0 

This equation is easily derived by considering the very 
first step in rapidity and with the help of the probabil- 
ity distribution (JTJ) . This is actually a version of what 
we call in QCD the Balitsky-Kovchegov (BK) equation 
[1, H3, Hl| associated to this simplified model. The phys- 
ical 5-matrix element in this framework, Sbk(jm)i is a 

particular value of Z, obtained from the above equation 

2 

by setting Xk^n = I and Xi — e~ a . The nonlinearity 
in Eq. ((4]) makes sure that S be unitary. We will define 
the scattering amplitudes properly only in Sec. IIII1 but 
before, we may summarize the main known properties of 
Sbk seen as a solution of Eq. (TJ| . 
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B. Basic properties of Sbk 

We know that Sbk nas the form of a wave front that 
travels towards larger values of \i\ under rapidity evo- 
lution. Its known properties are very well documented, 
and we refer the reader to the original papers in QCD 
[22l [23l [24| (see also Ref. [I^] for the pioneering calcu- 
lations, but without the connection to traveling waves) 
or to reviews in mathematical physics [2|| for the de- 
tails. Traveling waves have a phenomenological signature 
in high-energy scattering, called "geometric scaling" and 
found in deep- inelastic scattering data [13] • We give here 
the technical features of our particular model without 
much justification. 

When i and y are large, 



Sbk(2/,«) 



1 



-7o(i-I v ) 



(5) 



where X y is the position of the wave front at rapidity y. 
In the context of QCD, X y would be the logarithm of the 
squared saturation scale. 

As well-known, 70 and X y are determined from an ana- 
lysis of the linearized equation We insert Eq. ([5]) into 
Eq. ((4]) , and after linearization and some easy algebra, we 
arrive at the expression for the velocity V of a front of 
the form with decay rate 7: 



where 



X(7) = In 



V 



(6) 



1 - ei- 1 1 - e-T" 1 



is the characteristic function of the evolution kernel in 
Eq. ((4]). The relevant value of 7 at large rapidity is the 
one that minimizes V(j). This minimization gives 



70 = 0.607187- 



V Q = 1.02935 ■ 



(8) 



which are the decay rate and the velocity of the front at 
infinite rapidity. The asymptotics is approached as 1 



V(y) 



3 1 3 
Vb - ^ — + ^2 



2tt 



270 y 27 2 V X "( 7o ) y3/2 



(9) 



where the dots stand for subleading terms whose analyt- 
ical expression is not yet known. 



Starting from an initial condition of the form 

S(y = Q,i) = l + (e- a2 -l)S a, (10) 
that is to say, 

P(y = 0, n k ^ = and n k=0 = 1) = 1, 



P(y = 0, all other config.) = 0, 



(11) 



a front is formed (i.e. blackness is reached; S(y, 0) <C 1) 
after a rapidity of the order of 



Vf 



X (0) ln «2' 



(12) 



and it relaxes to its asymptotic shape up to a resolution 2 
a 2 after an evolution over 



1 2 



1 



27oV(7o) « 2 



(13) 



units of rapidity. 



III. UNITARY SCATTERING AMPLITUDES 

A. Formulation of scattering 

We consider the Fock state of a particle initially at 
position after a rapidity evolution y. Through the evo- 
lution, one eventually gets a system of {rij} particles. 
The probability that this system does not interact with 
a target consisting in a single particle at position i simply 
reads 

e~ a2n \ (14) 

This, of course, is also what we shall call the S-matrix 
element for forward elastic scattering. With a target that 
consists in a set of {rrij} particles on the sites indexed by 
j, the probability that there be no interaction reads 



n« 



(15) 



We have just assumed the complete independence of the 
individual scatterings, and that each of them occurs with 
probability a 2 if the objects in presence have the same 
position on the lattice, and with probability if this is 
not the case. Note that there is here a difference with the 



1 The first term Vq was already discussed by Gribov, Levin and 
Ryskin in QCD in the 80's y|. The second term was around 
in statistical physics since some time, see Ref. |2(I . but was re- 
discovered independently by Mueller and Triantafyllopoulos in 
QCD |25|. The third term was computed more recently [26l |. and 
adapted to QCD in Ref. [SI . 



2 We mean that the front is in its asymptotic shape (0 for S < 
1 — a 2 . Indeed, the asymptotic shape sets in first in the vicinity of 
the black region, and subsequently diffuses upwards. If one is not 
able to resolve details of size greater than a 2 , then it is enough 
that S look asymptotic in the above region. This distinction is 
relevant when one goes beyond the BK equation and one takes 
into account the fluctuations, and yn is then a physical relaxation 
"time" . 
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QCD dipole model since there, the interaction consists in 
at most one gluon exchange between each pair of dipoles, 
whereas with our choice any number of exchanges may 
occur. 

Let us consider two particles, initially at respective 
sites and i, which evolve into systems of nAy) and 
rrijii/) particles respectively on site j after a boost at ra- 
pidity y. In the Mueller approach [l3|, Qjj], the S-matrix 
for the scattering of these objects is defined, for low ra- 
pidities, as 

S(y,i) = iJjIe-^^N^di-'wj ( 16 ) 

where the average is taken over the realizations of the two 
systems. The i index is implicit in the r.h.s.: It is related 
to the initial condition that leads, after evolution, to the 
system {rrij}. a ranges between and 1 and defines the 
share of the rapidity between the two systems. When a — 
(or 1: the target and the projectile may be exchanged, 
hence there is a symmetry a — > 1 — a after the average 
over the realizations), the scattering occurs in the lab 
frame. In this case, the above formula gives the solution 
to the Balitsky-Kovchegov equation if the states evolve 
linearly: S — Sbk, where Sbk was discussed in Sec. Ill Bl 
When a — | instead, it means that the scattering takes 
place in the center-of-mass frame. 

By definition, S is unitary, whether or not there is a 
limit on the growth of the number of partons that inter- 
act. So a priori there is no violation of unitarity, and 
saturation is not necessary to prevent S from becoming 
negative. However, we know that the discussion is more 
subtle, and that the S-matrix defined like that violates 
boost invariance Q. Nevertheless, for small enough val- 
ues of y, there should be no need to put saturation effects 
in the evolution of the states. 

Let us discuss quantitatively the expected range in 
which saturation effects may be left out. In these consid- 
erations, we follow the discussion given by Mueller [l3|. 
To simplify, let us first go to the center-of-mass frame, 
that is, we set a = h. Imagine that we start from a 
situation in which both objects are at rest. When the 
two systems are gradually boosted in opposite directions 
in order to increase the center-of-mass energy of their 
eventual scattering, the probability of a high number of 
partons in each of them increases. At some point, when 
the unitarity limit is reached (as soon as the total ra- 
pidity is larger than yp), the probability that there be 
two or more interactions between them becomes of order 
1, while each of the Fock states remains relatively dilute 
(they contain only typically l/a particles each), in such 
a way that one can still consider their evolution as linear. 
Actually, because the systems share half of the rapidity, 
they become subject to nonlinear effects (in the form of 
internal rescatterings, for instance, or recombinations) at 
rapidity 2xyp- So this sets the limit beyond which multi- 
ple interactions between linearly-evolving Fock states are 
no longer sufficient to fully describe the scattering. The 



extension to an arbitrary value of a is straightforward, 
and we get the limit 



max(er, 1 — a) 

The center-of-mass frame (a = h) is naturally the most 
favorable case, with the highest limit on y, while in the 
lab frame, the wavefunction evolution can be linear only 
until blackness of the scattering amplitude is reached. 



B. Saturation 

So far, we have considered Fock-state evolution 
through splittings only (Sec. [TTJ) . A unitary S'-matrix 
was constructed from multiple scatterings in the center- 
of-mass frame with linearly-evolving states (Sec. IIII A[> . 
However, we know that theoretically, this procedure has 
a limited validity, that we have estimated as y < 2yp . In 
order to compare the results for scattering obtained in 
that framework to what would be obtained in the case 
of saturated states, we need to introduce a saturation 
mechanism. 

We can imagine different ways of enforcing saturation. 
The simplest way consists in forbidding new splittings 
to a site on which the number of particles has reached 
the value l/a 2 . We will call this method "veto". The 
resulting model is close to other models that have been 
studied before, e.g. in Ref. [28|. 

However, in order to get an approximately boost- 
invariant S'-matrix 3 , one needs instead to replace the 
splitting probability ^ at site i by ^e _Q ni , in such a way 
that for small rii compared to l/a 2 , the splitting prob- 
ability tends to the free splitting rate i, and for large 
7^, splittings are frozen. This prescription was proposed 
by Mueller and Salam [19( long ago and revisited more 
recently (29| . We will check numerically the approximate 
boost invariance. This method will be called "satura- 
tion". It is this one that makes sense in a field-theory 
framework where of course observables have to be inde- 
pendent of the frame in which the scattering is observed. 
Note that with this choice for the saturation mechanism, 
the resulting model is close to the one that was exten- 
sively studied in Ref. [29j |. except for the discretized ra- 
pidity. 

To our knowledge, the uniqueness of the prescription 
that leads to boost-invariant saturation has not been for- 
mally proven in models with spatial dimensions such as 
the one that we are building. It is also not known whether 
the prescription could be slightly modified in a way that 
approximately preserves boost invariance, the character- 



Since rapidity is discrete in this model, we can only have approx- 
imate boost invariance. We will observe the consequences of this 
incomplete realization of the symmetry in our numerical results. 
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istics of the amplitude being at the same time signifi- 
cantly changed. 

Finally, we note that the consequences of boost invari- 
ance (or rather of "projectile-target duality") has been 
thoroughly investigated very recently, also at a quite for- 
mal and rigorous level [3(J [3l| : The transformation that 
corresponds to exchanging the scattering particles was 
translated into an operation on the effective action of 
the corresponding model. So we could check the (ap- 
proximate) boost invariance of our model with the tech- 
nology developed in Refs. [2(| [Hj]. (We have not done 
so: We shall check boost invariance later on, but only 
numerically) . 



C. Expected properties of the S-matrix 

If the partonic evolution is linear (i.e. without satu- 
ration effects), and if the scattering takes place in the 
lab frame, then S = Sbk, and S exhibits the properties 
listed in Sec. HlBl 

We also know that, whatever the saturation mecha- 
nism is, as soon as the growth of the number of partons 
on each site is limited once it approaches 1 /a 2 , then the 
traveling wave solution for S is modified H2> HI, H3 • 
Essentially, its asymptotic velocity is corrected as follows 



V = V 



^ 2 "fox" ho) 



2(ln(l/a 2 ) + 31nln(l/a 2 )) 2 



(18) 



This equation contains the first term of an asymptotic 
expansion when ln(l/a 2 ) ^ 1 (It is actually valid up 
to 0(lnln(l/o; 2 )/ln 3 (l/Q; 2 )). We do not intend to go to 
extremely small values of a, thus the exact value of these 
asymptotics will not be probed. What is interesting and 
universal however is that V is less than V$. 

We will be able to extract more information from our 
Monte-Carlo, that we can compare with non-trivial and 
characteristic predictions made for reaction-diffusion sys- 
tems. In particular, since the whole scattering process is 
stochastic, there is a dispersion in the position of the 
front between different events. It reads [33| 



a 2 ) -a) 2 - t^lM 

{y) {v) ~31n 3 (l/a 2 ) 



y 



(19) 



While again the exact value will not be probed for it is too 
much asymptotic, the characteristic feature of the vari- 
ance of the position of the front is that it scales linearly 
with y. We recall that this linearly growing variance is 
at the origin of the phenomenon of "diffusive scaling" in 
the observable amplitude, which breaks geometric scal- 
ing predicted that comes out of the solution to the BK 
equation. 

As for the case when evolution does not include a sat- 
uration mechanism in the Fock states, it is difficult to 
figure out a priori what is going to happen. We will find 
out in the next section by performing numerical simula- 
tions of the different variants of the model. 



IV. NUMERICAL STUDY OF THE TOY 
MODEL 

We have implemented the model defined in Sees. HT1 and 
IIIII in the form of a Monte-Carlo event generator. There 
was no major difficulty to be overcome: The techniques 
that we have used are standard and the code can easily 
be reproduced by the reader. 

The evolution starts with one particle on a given site in 
each of the two systems. We evolve one or the other sys- 
tems by steps of one unit in y . We measure the position of 
the traveling wave front by searching, at each rapidity y, 
the rightmost site if for which S(y,ip) < 0.5. We then 
define the front position 1 y from a linear interpolation 
between ip and ip + 1. 

Let us move on to the results. 



A. S-matrix 

We compute the S'-matrix in 3 different variants of the 
model: bare multiple scatterings in the lab frame (which 
is the BK assumption; a = 0), the same but in the center- 
of-mass frame (which is Mueller's unitarization proce- 
dure; a = o) and multiple scatterings off boost-invariant 
saturated wave functions in the center-of-mass frame {S 
should be boost invariant in this case: We will check it 
later). 

Tuning the frame in the Monte-Carlo is technically 
easy: It is enough to share the rapidity evolution of the 
projectile and of the target proportionally to er. 

The results are presented in Fig. [2J We see that ra- 
pidity evolution starts with the blackening of the central 
region around i ~ 0. At y ~ 15 — 20 two traveling waves 
form, and propagate symmetrically towards larger (resp. 
smaller) values of i. In the following, our comments will 
always refer to the traveling wave that propagates along 
the positive i axis. 

We see that the fastest wave is observed in the BK 
model, while the saturation model gives rise to the slow- 
est one. We also observe that all waves get slanted during 
their propagation, except in the BK model. This is "dif- 
fusive scaling" , while the BK solutions exhibit geometric 
scaling. 

In the next paragraphs, we go deeper into the analysis 
of this calculation by focussing on the statistics of the 
position of the front: Its mean velocity and its variance, 
for the different variants of the model. 



B. BK equation 

We go to the lab frame, that is, we put all the rapidity 
evolution in one of the objects, while the other one is 
at rest. Multiple scatterings off a system that does not 
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FIG. 2: ^-matrix element as a function of the spatial coordinate for y — 5, 10, 20, 30 (from the center of the figure towards 
the outskirts), l/a 2 — 20. For low rapidities (y = 5 and y = 10), the evolution is linear (of BFKL type) and the number of 
particles grows exponentially. At y = 20, the unitarity limit has been reached in all calculations. Later (y — 30), the traveling 
waves are formed and propagate outwards. Geometric scaling violations are seen in the tilting of the fronts, except for BK-type 
scatterings. 



saturate gives a S-matrix that solves the BK equation, 4 
and hence that has the form of a traveling wave whose 
characteristics are given in Sec. Ill Bl The velocity of the 
front is presented in Fig. [3] for l/a 2 = 20. In order to 
interpret the results, it is useful to compute the numerical 
values of y F and yn in Eq. (fT2"]) and (flU]) : 

y F ~ 7.4, y R ~ 2.2 for l/a 2 = 20. (20) 

We see in Fig. [3] that first the velocity is 0, which cor- 
responds to the phase in which the parton numbers are 
building up through a linear evolution. (The linear phase 
is named in QCD after Balitsky, Fadin, Kuraev and Lipa- 
tov (BFKL) [37|). Later, after about 10 units of rapidity 
(which is the order of magnitude of y F ), a sharp peak 
appears and decays over a few units of y. This peak may 
be interpreted as follows: In the initial stage of the evo- 
lution, when the front is building up, its shape is steeper 
than the asymptotic shape (O in the region where the 
position is measured (S ~ 0.5), and hence its velocity is 
larger than Vq . Then, the front relaxes to its asymptotic 
shape in that region, which takes of the order of yn steps 
of rapidity. In the final stage, for y 3> 10, the asymp- 
totic velocity is approached in an algebraic way, in good 



4 Several groups have solved the exact BK equation numerically 
in QCD, see for example [2Stl35l. [36l]. A good code (BKsolver) is 
publicly available at http://www.isv.uu.se/~enberg/BK/ 



agreement with the theoretical expectations of Eq. ([9]). 



C. Front velocity in different models 

We want to compare bare unitarization through mul- 
tiple scatterings in different frames to the variant of the 
model in which saturation is included in a boost-invariant 
way. The front velocity for different schemes and l/a 2 
set to 20 is shown in Fig. Q] 

Bare multiple scatterings lead to a dip around y ~ 20 
whose depth is maximum for a = i. At large y, the ve- 
locity Vq is reached algebraically, whatever the frame is. 
At any finite y, this model gives rise to different front so- 
lutions in different frames: Boost invariance is manifestly 
badly broken. 

Solutions with saturation in the evolution look like ex- 
pected: The velocity rapidly reaches a plateau, with a 
value that is lower than that of the asymptotic BK ve- 
locity. Theoretically, this should happen after typically 
Vf + Mr units of rapidity. For low values of l/a 2 , the 
formula giving yn is expected to need large corrections, 
and so the numerical value of yn is not trustable. The 
value of the asymptotic velocity is equal in all frames, 
and except for the lab frame, all curves superimpose rea- 
sonably well during the whole rapidity evolution. We 
deduce that boost invariance is quite well preserved with 
our saturation solution, although not perfectly. We at- 
tribute the lack of complete superposition of the curves 
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FIG. 3: Front velocity for the solution of the equivalent BK 
equation as a function of the rapidity, compared to the theo- 
retical expectations: the different theoretical curves represent 
1,2 or 3 terms of Eq. ©. 1/a 2 is set to 20. 

to the explicit breaking of boost-invariance due to our 
discretization in rapidity. 

The same comments are true for larger values of 1/a 2 , 
see Fig. [5] and [6] It is useful to compute yp and ya also 
in that case: 

y F ~ 18.7, y B ~ 17.0 for 1/a 2 = 2 x 10 3 , 

y y ' ' (21) 

y F ~ 30.1, y R ~ 43.9 for 1/a 2 = 2 x 10°. 

We check that a non-zero velocity appears for y > yp. 
Actually, yF given in (|2ip systematically underestimates 
the actual value of the rapidity at which a front is formed. 

So far, we see in Fig. I4I5I6I that the characteristics of 
the front agree in the different models formulated in the 
center-of-mass frame only in a very small interval of ra- 
pidity, roughly consistent with the theoretical estimate 
y < 2yp. Except for large values of 1/a 2 , the front ve- 
locity in the boost-invariant scheme is not reached by 
bare multiple scatterings in the center-of-mass frame. 

It is instructive to also study the variance of the po- 
sition of the front in the different variants of the model: 
We analyze this quantity in the next paragraph. 

D. Variance of the front position 

We now turn to commenting on the variance of the 
position of the front. This is shown in Figs. I7I8I9I for 
different values of 1/a 2 . 

First, we notice a sharp difference between saturation 
and bare multiple scatterings. The saturation scheme 



FIG. 4: Front velocity for different unitarization schemes: 
Multiple scatterings (upper bunch of curves), and boost- 
invariant saturation (lower curves). Four distinct frames 
are considered: a = 0.5,0.25,0.1,0, respectively denoted by 
"COM", "1:3", "1:9" and "LAB". 1/a 2 is set to 20 in this 
figure. 



leads to linearly increasing fluctuations as soon as the 
front is formed, in qualitative agreement with Eq. p9[) . 
By contrast, the bare multiple-scattering schemes lead to 
fluctuations that slow down with y. However, there is a 
range in rapidity in which the variances agree in the dif- 
ferent models and in all frames, except for the lab frame 
which has systematically less fluctuations. The agree- 
ment is particularly striking at large 1/a 2 . The range in 
which the models match can be estimated as y < 2yp for 
the center-of-mass frame, if yp is taken to the rapidity 
at which the front is effectively formed (around the lo- 
cal maximum in the variance) rather than the numerical 
estimate done before. 

The fact that the calculations in the lab frame lead to 
different front velocities and quite different fluctuations 
should maybe not come as a surprise, since the way in 
which we treat this frame is very special: Indeed, we do 
not allow at all fluctuations in one of the incoming ob- 
jects, which is quite unphysical (quantum objects should 
be allowed to fluctuate even if they have a vanishing ra- 
pidity), and which is likely to reduce the event-by-event 
fluctuations that are observed in the scattering of the 
objects. 

As was recalled in Sec. IIII CI the linear growth of the 
variance of the front position is characteristic of reaction- 
diffusion processes, and is very well seen in the model 
with saturation, already for moderately small rapidities. 
The fact that this linear behavior is reproduced by multi- 
ple scatterings in the center-of-mass frame (see especially 
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FIG. 5: The same as in Fig. H but for 1/a 2 = 2000. 



1.15 



1.05 



0.0". 



0.9 



U.N.", 



-COM 
-LAB 





FIG. 7: Dispersion of the position of the front. The curves 
that converge to a straight line at large y correspond to Fock 
states evolved with a kernel that includes boost-invariant sat- 
uration, while the ones that flatten correspond to multiple- 
scattering unitarization. 1/a 2 is set to 20. Inset: Zoom into 
the region of low rapidity for the calculations in the center- 
of-mass frame. 



the models agree in the case l/a 2 = 2x 10 5 , the asymp- 
totic slope of the variance has not been reached yet. So 
except for a model in which ur ^ uf and yp ^> 1, which 
are two conditions that are difficult to realize in actual 
models, it is very difficult to make convincing statements 
on the ability of bare multiple scatterings to mimic a 
reaction-diffusion process. 

So far, we have considered boost-invariant saturation 
only. Boost invariance is a basic requirement for the 
model to be consistent with field theory. We wish how- 
ever to compare to a non boost-invariant scheme. 



20 40 60 80 100 120 140 

y 

FIG. 6: The same as in Fig. H but for 1/a 2 = 2 x 10 5 and 
only two frames (center-of-mass and lab) are used. 



Fig. [8]) over some finite range of validity may be some ev- 
idence in favor of the reaction-diffusion interpretation of 
high energy scattering. But admittedly, this linear be- 
havior is seen in a very small range for small 1/a 2 . On 
the other hand, for large 1/a 2 , yn > yF and thus the 
front has not properly relaxed before genuine saturation 
effects should be taken taken into account, for y ~ 2yp- 
As a matter of fact, we see that in the range in which 



E. Other (non boost-invariant) saturation scheme 

We adopt the alternative scheme of saturation which 
consists in vetoing further particle splittings to a site as 
soon as the number of particles on this very site has 
reached the value 1/a 2 (see Sec. IIII B|) . The small-a 
asymptotics of the statistics of the front position should 
be insensitive to the exact way how saturation occurs. 
However, subleading effects at finite a have no reason to 
be identical. 

We see indeed in Fig. [TU] that the asymptotic velocity 
is higher for this scheme than for boost-invariant satu- 
ration. Moreover, multiple-scattering unitarization leads 
to a velocity that, at low y, is closer to the one obtained 
from the veto procedure. 
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FIG. 8: The same as in Fig. [7] but for 1/a 2 = 2000. 
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FIG. 9: The same as in Fig. [7] but for 1/a 2 = 2 x 10 5 



However, if one observes the variance of the position 
of the front (Fig. [TTj) , one sees that the large- y slope is 
quite different between the veto scheme and the boost- 
invariant scheme. Around y ~ 50, it is clear that the 
unitarization scheme is much closer to boost-invariant 
saturation. On the other hand, the large-?/ slope of the 
variance in the boost-invariant scheme looks like the con- 
tinuation of the slope in the domain in which they agree. 

In this sense, boost-invariant saturation is the natu- 
ral continuation of the multiple-scattering unitarization 



FIG. 10: Comparison of the front velocities obtained with a 
linear evolution and multiple-scatterings in the COM system, 
and with the two different ways of adding parton saturation 
(veto and boost-invariant saturation). 1/a 2 = 2000. 



procedure, since it "knows" about the fluctuations of sat- 
urated Fock states at large values of the rapidity. The 
whole difficulty would be to find how, in practice, to per- 
form this continuation. 



V. ALTERNATIVE CALCULATION: 
RESUMMATION OF ASYMPTOTIC SERIES 

A. Pomeron-loop expansion 

In the framework of the approximations of high-energy 
scattering that we are considering (multiple scatterings of 
unsaturated Fock states in the center-of-mass frame) , the 
scattering cross section may be computed from Eq. (|16p . 
That expression may also be further expanded in powers 
of a 2 . We get the following series: 



2\k 



fc=0 



k\ 



/2)m i (y/2) 



(22) 



k is the number of Pomerons exchanged between the tar- 
get and the projectile: k = 1 is the tree-level (BFKL) 
term, k = 2 is a one-loop contribution, and so on. 

The series that has been obtained is actually a diver- 
gent series: The term of order k behaves like kl because 
essentially, (n k ) ~ fc!. The problem is now the following: 
Assuming that we know the first terms in this divergent 
series, can we get an estimate of the fully resummed se- 
ries? 



10 



20 - 



Saturation 
Veto 

Multiple scatterings 




FIG. 11: The same as in Fig, 1101 but now the variance of the 
position of the front is shown instead of the velocity. Inset: 
zoom into the region in which the three calculations almost 
agree. 



B. Overview of the resummation method 

We have to perform an infinite sum, where only a fixed 
number of terms is known. There arc different possibili- 
ties to estimate the limit of this sum based on these re- 
stricted pieces of information. Depending on the asymp- 
totic behavior of the partial sums, literally summing the 
first elements of a series in many cases is too slow or even 
does not converge, as in our case. In physics, the best 
known techniques to deal with divergent series are Borcl 
summation [381 ] and Pade approximation [391 ] . 

On the other hand, there exists nowadays an enormous 
number of non-linear sequence transformations which can 
accelerate the convergence of convergent series and which 
have also shown to be very efficient in summing divergent 
series. Particularly suited is a class of such transforma- 
tions introduced by Levin [4(| ■ For an introduction to the 
topic we refer the reader to Ref. [4l|, \^ and references 
therein. Here we just sketch the main formulae. 

Consider the partial sums s n with limit (or "antilimit" , 
which is how the resummed value of a formally divergent 
sum is called) s and remainder R n : 



z=0 



(23) 



The aim now is to find a new sequence of partial sums s' n 
such that R n /R n — > for n — > oo. An important building 
block of a specific sequence transformation is a remainder 
estimate Lo n which should reflect the behavior of the exact 
remainder. Since the remainder estimate only describes 



the leading behavior of the exact remainder, we write 



Rn /^n^n? 



(24) 



where ji n is of the order of 1 and converges to some 
(unknown) number. The second specific ingredient is 
a set of functions tjji(n) on which the /j, n are decom- 
posed. For the original Levin-transformation [4(| the set 
ij}i(n) = (n+[3)~ l was chosen, where (3 is an arbitrary real 
positive parameter which in general is set to 1 (We will 
however choose a different value of (3 in our application 
of the method). We write 



oo 

y^c^(n) 



oo. 



(25) 



Of course, the coefficients are unknown, but if we trun- 
cate the sum in Eq. (125[) . we can interpret Eq. (|2"3"]1 as a 
model sequence 



fc-i 



(26) 



i=0 



Inserting the values of the partial sums at hand for er m , 
we have for m S {0, 1, . . . , k} a system of linear equa- 
tions which can be solved exactly for a by Cramer's rule. 
By a recursive approach one can deduce a compact ex- 
pression for a which circumvents the evaluation of large 
determinants: 



-i , {k~ 

=0 A 0,i 



-i 

=0 A 0,i 



(27) 



where the coefficients \ have to be calculated for a 
given set of functions ipi(n) and are independent of the 
concrete remainder estimate. At this general level, rigor- 
ous mathematical statements about the convergence are 
hardly possible, but for well-posed expansions as they ap- 
j>ear in physical problems a converges to s when k — > oo 

For the remainder estimate, Levin introduced three 
variants (see Tab. |T|. The t- variant is adapted for the 
case of linear convergence, while the u- and u-variants 
shall be also usable for logarithmic convergence. The d- 
variant has been proposed in Ref. [44| especially for alter- 
nating logarithmically convergent series. The c-variant 
has been designed for the same purpose [i^] • All of them 
have been used for the summation of divergent sums as 
well. In face of these preliminary considerations, there is 
no strict rule which remainder estimate one has to use. 
For our problem, it turned out that variants c and d are 
less successful to the other ones, where u and v give the 
most stable results. 

The convergence significantly improves if one compiles 
inverse factorials (or more accurately Pochhammer sym- 
bols) instead of inverse powers to the function set il>i (n) 
leading to the Weniger S- or A^-transformations [43j |. 
The explicit expressions for the functions ipi(n) and the 
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(n + /3)a„ 


a n a n4-l 




(inan + 2 







TABLE I: Remainder estimates for different sequence trans- 
formations. 

coefficients used in Eq. ([2"7| are given in Tab. [Til 
They can be combined with the same remainder esti- 
mates and are then labeled (analog to Levin-transforms, 
see Tab. UJ r, y, <fi, 5, \ f° r the iS-transformations, and 
T, Y, <I>, A, X for the .M-transformations. While the 
.M-transformations are of no use for us, as long as we 
stay with (3=1, the iS-transformations provide good re- 
sults with the same ranking concerning the remainder 
estimates, i.e. the most stable results are obtained us- 
ing the (^-transformation. But the remainder estimate is 
not pure trial-and-error. From the knowledge about our 
sum it is clear that the remainder estimate t correctly 
describes the behavior of the sum. 

What goes wrong? Usually these sequence transfor- 
mations are used when the available elements are known 
with high precision. By contrast, we calculate our ele- 
ments by a Monte Carlo simulation in which the higher 
moments are afflicted by larger relative statistical error 
than the smaller ones. Moreover, these higher moments 
cause a statistical error by their mere size compared to 
the other moments when implemented on a computer 
with limited precision. The sequence transformations 
discussed so far emphasize these larger moments assum- 
ing that they are already closer to the limit. If we re- 
frain from a too large impact of these larger moments, 
we can shift to a more balanced combination of the mo- 
ments by increasing (3 in the coefficients X^\. In the 
limit (3 — > oo one would obtain a symmetric weighting of 
large and small moments X^] = (— known as the 
Drummond-transformation [43L |45| . 

type Mn) A^j 

L-in ^ (-iy(*)£±£±^ 

WenigerS ^ (-l)'O ^ 

WenigerM (-!)«(*) {l^Zg^ 

TABLE II: A-coefficients for different sequence transforma- 
tions. 

Finally we use the remainder estimate of the Weniger 
r- variant with f3 = 100. 

Beside these Levin-type transformation there exist also 
many other schemes which are less generally applicable. 
They cannot be written in the form of Eq. (|27|l but are 
given by a recursive definition. From these we also tried 
the e algorithm [46[ and thereby Shanks transformation 
E3, the Aitken A 2 process in its iterated form (43l . 
[49L the p J transformation [5(|, Brezinski's 9 algorithm 
[5ll | and its iterated formulation (43j . 



C. Numerical implementation and results 

We wish to apply the resummation methods described 
above to the computation of S in Eq. ([22]) . 

We need to compute numerically the first few terms 
of the series ([22)) . These are actually proportional to the 
moments of the one-Pomeron exchange amplitude, which 
reads 

a 2 J2^(v/ 2 )^(y/ 2 l (28) 
j 

where {rij} and {rrij} are realizations of systems of par- 
ticles (linearly) evolved by the Monte Carlo algorithm 
described above. 

The numerical implementation of this calculation is 
straightforward. For each event, the one-Pomeron am- 
plitude is evaluated for all values of i and y. Then, its 
fc-th power is computed, and the average over events is 
performed. We then apply the resummation methods 
described above to get the result for S. 

While for the calculations of Sec. IIVI one million of 
events were enough, here, we need hundreds of millions 
of them. Indeed, it turns out that statistical errors are 
large for higher moments of Eq. ([28]) . 

In practice, we generate 5 x 10 8 events. We consider 
the resummation of k = 5, 6 • • ■ terms, up to 30, for dif- 
ferent values of a 2 . We consider that this is the present 
technical frontier, since a few months of calculation were 
already needed to achieve this number of events. 

The result of the resummation for S is shown in Fig. [T^l 
for 3 different rapidities and for different number of terms 
k, as a function of the site index. First of all, we see a 
good convergence of S when the number of terms that are 
taken into account increases. Second, we see that for y = 
15 and y = 25, S obtained from this resummation looks 
exactly like the one obtained from multiple scatterings in 
the center-of-mass frame. We see that for y = 50 instead, 
the result of the resummation does not coincide with any 
of the previous calculations. 

In order to have a more synthetic estimate of the qual- 
ity of the resummation, we can compute the velocity of 
the front and compare it to the calculations in Sec. IIVI 
We see in Fig. [T3] that all calculations match for low y. 
When more terms are taken into account (larger value of 
k), the domain in which the statistical and asymptotic 
series calculations agree extends towards larger values of 
k. For k = 30, the agreement is very good up to values 
of y of the order of 25. 

To gain confidence in the stability of our resumma- 
tion, we can compare two out of the many resummation 
methods described in the previous section in Fig. [Ml We 
see a very good agreement for all k. We conclude that 
the discrepancy at large k with the calculation of Sec. ITVl 
is not due to a failure of the resummation method, but 
rather to a lack of the relevant information, which would 
be contained in higher-order terms k > 30. We would 
also like to show a method which does not converge very 
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FIG. 12: Resummed S'-matrix for 3 rapidities. Different calculations are compared: multiple scatterings in the lab and COM 
frames and boost-invariant saturation (points). The three curves correspond to the resummation of the asymptotic series using 
the Weniger ^-transform for 5,15 and 25 terms respectively. 




FIG. 13: Front velocity as a function of the total rapidity from 
the resummation of the asymptotic series (bunch of curves in 
continuous line: the highest up corresponds to k = 5, the 
lowest to k = 30) compared to the different calculations. 
l/a 2 = 20. 



FIG. 14: Comparison of two different resummation methods 
for the front velocity: The Weniger 5-transform (r version; 
method I), and the Weniger .M-transform (T version; method 
II). In each case, the parameter f3 is 100. Low fc's lead to 
higher values of the velocity in the region y ~ 20. 



well: Therefore, we do the same but setting the parame- for high values of k does not give a meaningful result, 
ter (3 to 10 (instead of the higher value 100 that we have More statistics would however help (more events gener- 
chosen so far). We see in Fig. [L31 that the resummation ated, that is, a better accuracy in all terms). For yet 
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FIG. 15: Resummation using a lower value of the parameter 
f3, namely ft = 10. 
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FIG. 16: The same as in Fig. [T3 but with 1/a 2 = 2000. 



lower values of /3, such as (3 = 1, the result would even 
be worse. 

Finally, we set 1/a 2 = 2000, and we see again the 
same features, and in particular a good agreement with 
the calculations of Sec. IIV1 this time up to y ~ 35 (see 
Fig.[H|). 

Note that we have limited our study to the front ve- 
locity averaged over events, in the statistical language of 



Sec. IIII1 The variance of the position of the front would 
require a special calculation, that we have not done. 



VI. CONCLUSION AND OUTLOOK 

The goal of this paper was twofold. First, to test the 
original unitarization procedure through multiple scat- 
terings, without explicit saturation, proposed by Mueller 
in the context of the color dipole model, in the light of the 
new understanding of unitarization gained from the anal- 
ogy with reaction-diffusion systems. Second, to try and 
reproduce these results from an expansion in the number 
of Pomerons that are exchanged. 

We have shown empirically that the original formula- 
tion of unitarization by Mueller is successful within the 
(limited) range of validity that had been assigned to it. 
The traveling waves exhibit a front velocity that is infe- 
rior to the one that would be expected for a system with 
no saturation mechanism at all, and the event-by-event 
dispersion in the position of the front grows linearly with 
the rapidity, as expected for reaction-diffusion systems. 

We have resummed the multiple scatterings in two 
ways. The one that was used for example by Salam in 
his Monte-Carlo, which consists in averaging over events 
the scattering matrix, and another one which relies on 
an expansion of the S'-matrix in powers of a 2 . A pri- 
ori, it was not completely obvious that these two ways of 
computing the unitarized dipole-dipole scattering cross 
section would lead to the same result for the S'-matrix, 
since one sums the defining series in a very different or- 
der in both cases. But the fact that we eventually get 
the same answer is reasonable. 

The resummation tools that we have set up and tested 
in the simple toy model studied here will be useful in 
the future. Indeed, if we were able to compute order by 
order the unitarity corrections, either in this toy model 
or in full QCD [12, H3], we would be confronted to the 
task of summing the resulting series, which has the struc- 
ture of the one that we have studied. In general, there 
is no reason why there should be a simple formulation 
like Eq. (|TB|) for scattering amplitudes, and field theory 
would naturally lead to an asymptotic series in powers of 
a 2 . We see that for moderately small rapidities and a, 
resumming the asymptotic series numerically is doable, 
although very difficult due to the large number of events 
that one has to generate in order to achieve a sufficient 
accuracy. 

In our next publication, we intend to systematically 
study the corrections to the Mueller formulation of uni- 
tarization, in the framework of our toy model. We are 
curious to find out whether boost-invariant saturation 
may be obtained through splittings in the Fock states 
and scatterings only, that is, without any explicit refer- 
ence to a saturation mechanism in the formulation of the 
model. This is of course a crucial question for QCD, since 
saturation seems so difficult to formulate in a practical 
way. 
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Another prospect would be to find an analytical ex- 
pression for S within the Mueller formulation of unita- 
rization: Since the result contains information on satura- 
tion, an analytical expression would help. But this would 
require the knowledge of all moments of the particle num- 
bers, which are given by the solution of the equivalent of 
the Balitsky hierarchy. Even within simple models, such 
an achievement does not seem to be at hand yet. 
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